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Abstract 



I previously used Burgers' equation to introduce a new method 
of numerical discretisation of pdes. The analysis is based upon cen- 
tre manifold theory so we are assured that the discretisation accu- 
rately models all the processes and their subgrid scale interactions. 
Here I show how boundaries to the physical domain may be natu- 
rally incorporated into the numerical modelling of Burgers' equation. 
We investigate Neumann and Dirichlet boundary conditions. As well 
as modelling the nonlinear advection, the method naturally derives 
symmetric matrices with constant bandwidth to correspond to the 
self-adjoint diffusion operator. The techniques developed here may 
be used to accurately model the nonlinear evolution of quite general 
spatio-temporal dynamical systems on bounded domains. 
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1 Introduction 

We discuss the holistic discretisation of boundary conditions for partial dif- 
ferential equations. Holistic discretisation || is based upon the support cen- 
tre manifold theory gives to the nonlinear dynamics on finite grid spacing. 
We expect such discretisation will have good stability and high accuracy on 
coarse grids because it systematically accounts for subgrid scale interactions. 
To date we have considered periodic problems |6], [1], an d their initial con- 
ditions J7|. Here we show how to incorporate different boundary conditions 
into the analysis. 

As an illustrative example we restrict attention to the one-dimensional 
spatial discretisation of Burgers' equation 

du du d 2 u 
dt dx dx 2 ' 

which contains the important mechanisms of diffusion and nonlinear advec- 
tion. As example boundary conditions we consider the important cases of 
Dirichlet, §131, and Neumann boundary conditions, §^j. The same techniques 
may be easily extended to other partial differential equations, other boundary 
conditions and higher spatial dimensions. 

The centre manifold analysis is based upon dividing the domain into finite 
sized elements, each separated from their two neighbours by specially crafted 
artificial internal boundary conditions (iBCs). The form of the IBCs (|3j-§|) 
generates a discretisation in the interior of the domain (§|2|) which is not only 
linearly consistent, as proved in ||, but which appears to be nonlinearly con- 
sistent to high order. This observation is based upon this analysis of Burgers' 
equation and work in progress on the Kuramoto-Sivashinsky equation; fur- 
ther research is needed to prove it in general. 

But the focus here is on the discretisation near the boundary of the do- 
main. Boundary conditions are easily incorporated simply by replacing an 
IBC of an end element by a variant of the actual boundary condition: the 
Dirichlet boundary condition of fixed field u is implemented as (|9D; the Neu- 
mann boundary condition of fixed flux u x is (|TT[) . The computer algebra 
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of §|A| readily computes the effect these boundary conditions have on the 
discretisation near the boundary. Novel features of this analysis are the fol- 
lowing: 

• theoretical support is provided on finite grid spacing h as explained 
in IBB e.g.]; 

• the discretisation has a consistent bandwidth across the whole domain 
including near the boundaries; 

• high order discretisations are expressed in terms of the matrices of the 
basic, second order, centred difference operators, leading to appropri- 
ately symmetric discretisations of the diffusion operator; 

• time variations of boundary values, say a, not only has a direct effect 
but also involves time derivative factors, such as d, which are more im- 
portant on coarser grids due to the time it takes changes in a boundary 
value to diffuse into the element. 

2 Holistic discretisation of Burgers' equation 
in the interior 

We consider Burgers' equation ([!]) on some finite domain in x. Place grid 
points Xj equi-spaced across the domain with constant spacing h, and corre- 
spondingly define Uj = u(xj, t), that is, Uj is the evolving field u evaluated at 
the j'th grid point. Divide the domain into m elements with a grid point at 
the centre of each. Following H to apply centre manifold theory we distin- 
guish each element using internal boundary conditions (iBCs) in the discrete 
form 

IJ>xfixVj(x,t) = 7/i5 Uj and S^.Vj(x,t) = ■j5 2 Uj , (2) 

evaluated at x = Xj where Vj denotes the field in the jth element. As shown 
in ||, this particular choice of IBCs ensures that the resultant finite difference 
scheme is consistent to high order in h as the grid size h — > . In this section 
we repeat the construction of the holistic discretisation, following ||, away 
from the domain boundaries via centre manifold theory but with a new and 
convenient form of these discrete IBCs. 

In actually developing finite difference models the IBCs may take any of 
many equivalent forms || e.g.]. In later sections we investigate the discreti- 
sation near a boundary of the domain: the element adjacent a boundary has 
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one real boundary and one artificial internal boundary. Thus it is appro- 
priate to rewrite the two ibcs in (^) in the form of two conditions, one at 
the left edge of each internal element, and one at the right. Recall that the 
difference operators [i ± \8 — E^ 1 ! 2 [|3], p65, e.g.] so that to the first IBC 
in (H) add/subtract half the second to give the equivalent ibcs 

5 x Vj(x,t) = 7o"«f_i/2 at x — Xj-1/2 , (3) 
and 5 x Vj(x, t) = 7^+1/2 at x — Xj+1/2 ■ (4) 

The IBC (|4]) is to hold at the right-hand side of each element and the IBC (|3]) 
is to hold at the left. The introduced parameter 7, when non-zero, couples 
each element to its neighbour. When 7 = these ibcs effectively insulate 
each element from its neighbours and forms the basis of the centre manifold 
analysis; whereas when 7 = 1 they assert that the field in the j'th element 
when extended into the surrounding elements has the same differences cen- 
tred across each internal boundary as that given by the grid point values. 
Thus evaluating the model at 7 = 1 forms the relevant discretisation. These 
ibcs apply to all elements except for the leftmost and rightmost elements, 
the ones adjacent to the boundary: in the leftmost element the left-hand 
IBC, (||) with j = 1, is replaced by the actual boundary condition; in the 
rightmost element the right-hand IBC, (|j) with j = m, is replaced by the 
actual boundary condition. In this paper I analyse the left boundary of the 
domain — discretisations near the right boundary are similar by symmetry. 

In the interior of the domain the boundaries of the domain have no in- 
fluence upon the discretisation. Each element in the interior has ibcs (|3]) 
and (U). Executing the computer algebra program in Appendix |A], adapted 
from H, the subgrid structure of the solution field u(x,t) is, in terms of 
£ = (x-Xj)/h, 



Vj = Uj 



+ lhl(e-0u J S\ + O(\\uf + 1 3 ). (5) 

The corresponding evolution on this centre manifold forms the holistic dis- 
cretisation in the interior: 



1 

1 



x2 T jr4 1 7 3 c-6 
70 Uj Uj H Uj 

' J 12 3 90 3 



-~h U3 



^{ibuj — ^-fi5 3 Uj 



7 2 

+ ^Uj + 5% n5uj) 



+ ^S + <?(HI 4 + 7 4 ). (6) 
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As discussed previously || when the coupling parameter 7 = 1: the first 
line gives successive approximations to the diffusion term u xx ; the second 
line gives approximations to the nonlinear advection; and lastly, the cubic 
nonlinear term, u 2 S 2 Uj in the last line above, accounts for subgrid scale in- 
teractions between the advection and diffusion and acts to stabilise the nu- 
merical model. Higher-order terms are easily found by the computer algebra 
of Appendix [A] but for clarity are not presented here. In the limit as the grid 
spacing h — > higher-order discretisations have the equivalent pde 



du 
dt 



7 



du d 2 u 



-u- 



dx dx 2 

A' 
720 

+ 15u x u 



h 2 , 
+ ^7(1-7) [ 



U, 



2 5 WW x^ XXX I '^'£C3/ 



7<Jf|7(l 7) ( ^ U X U XX 9llU xx <U<JUlU, X U, XXX I lUU. n U m 

2u 2 u xxxx ) + (1 - 4 7 )(2^n - 6ud 5 x u + hu 2 d x u 



0(ll«l| 4 ,7 4 ,^ 



(7) 



See that upon substituting 7 = 1 to recover a discretisation for the Burgers' 
equation (|I|), we would find an equivalent PDE to an error 0(h 6 , ||w|| 4 ). Anal- 
ogous lower order discretisations are obtained, as promised by the analysis 
in 0] , when we truncate the discretisation (^) to lower orders in the coupling 
parameter 7. But observe the new feature that the discrete form (|||-§]) of 
the ibcs lead to discretisations which are not only linearly consistent, as as- 
sured in m, but are also nonlinearly consistent. Further research is needed 
to establish nonlinear consistency in general. 



3 Dirichlet boundary conditions applied at a 
grid point 

Consider the case of Dirichlet conditions on the boundary of the domain of 
prescribed u at a grid point; without loss of generality say 

u = a(t) at x = Xo . (8) 

This boundary condition is included in the analysis simply by replacing the 
left-hand ibc in the leftmost element, (|^) with j — 1, by (as if uq = a in (§)) 

S x Vi(x, t) = 7 (w 1 — a) at x = X\/ 2 ■ (9) 

(Implicitly the first element then extends from xo = X\ — h to x\ + h/2.) 
When the coupling parameter 7 = this ibc effectively insulates the first 
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element from the conditions at the domain boundary. However, when 7 = 1, 
since Vj(xj, t) = Uj, this reduces to (||) by requiring vi(x , t) = a . Hence, the 
centre manifold derivation is based on 7 = as explained in || and evaluated 
at 7 = 1 to obtain a discretisation of Burgers' equation. 

We then solve for the subgrid fields in the elements near the boundary, 
j = 1, 2, . . . . The computer algebra program in Appendix ^ implements this 
boundary condition when dirichlet :=1 . The influence of this specified 
boundary value affects a number of elements near the boundary equal to 
the order of 7 retained in the analysis, nothing else: in the interior the 
discretisation is (Eh; whereas for elements the near the boundary and for 



errors (9(||it|| 4 + 7 4 ) we find the evolution to be of the form 



II: l 

u 3 



fc 
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fc 
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(10) 



where U = diag(«i, U2, 1*3) is the diagonal matrix of grid velocities and the 
three parts of the right-hand side represent respectively the discretisation 
of the diffusion, the nonlinear advection, and the leading order interaction 
between advection and diffusion. These parts include the forcing due to the 
time dependent boundary value a. Here the various terms are found to be: 
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and B = j^D + 0(7). Denote by D the matrix appearing above in T> lin- 
ear in 7 and then row extended across the interior of the domain: D is the 
matrix of the second-order centred approximation to the second derivative, 
5 2 . Observe that the order 7 2 and 7 s matrices in T> simply correspond to O 2 
and D 3 . Thus the discretisation of the diffusion term, across the entire do- 
main including the near boundary elements, is simply 7<^ 2 — if^ 4 + go^ 6 seen in 
the first line of the interior discretisation, (^), but with the matrix D replac- 
ing the centred difference 5 2 . Thus this approach generates an appropriately 
symmetric discretisation of the self-adjoint diffusion term u xx . The nonlin- 
ear stabilisation term, U 2 Bu, also corresponds to replacing 5 2 in by D. 
Similarly, denote by C the matrix appearing above in C linear in 7 (including 
the factor |) and then row extended across the interior of the domain: C is 
the matrix of the second-order centred approximation to the first derivative, 
fi5. Observe that the order 7 2 matrix in C is (DC + CD)/2 corresponding to 
the average of different permutations of the centred difference operators /i<5 3 . 
Thus again the discretisation of the advection terms across the entire domain 
is the interior discretisation, 7/1$ — -^(fi55 2 + 5 2 p>5)/2, with the truncated C 
and D replacing the difference operators fiS and S 2 respectively. Additional 
terms represented by g c for the discretisation of the advection correspond to 
the nonlinear higher order term ^(fi 2 Uj fJ>5 3 Uj + 5 Uj fiSuj) in the model (^]). 
These are pleasing patterns. 

Also see that by truncating to a fixed power of the inter-element coupling 
parameter 7 we will obtain a discretisation that has constant bandwidth 
across the whole domain. This constant bandwidth will always be derived 
in this holistic approach, see for another example the Neumann boundary 
conditions in the next section, because the truncation at a fixed power of the 
inter-element coupling 7 controls how many neighbouring elements interact 
with any given element. Although there is a lower order (in h) of consistency 
near the domain boundaries, the support by centre manifold theory is the 
same across the whole domain. As discussed in |J, the reason for this support 
is that the theory applies to the solutions of Burgers' pde (|TJ) in the entire 
domain, not just in some locale. 

Now consider the forcing from the boundary. Using a = (a, h 2 a), it is 



Tf. _|_ ri 7 7 2 7 3 



6 ~ 18 12 45 112 



iL _ "if- I 7 2 I 7 3 

12 45 ' 90 ^ 140 



a + 0( 7 4 ,a), 



(14) 
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(15) 
(16) 



See that this holistic approach not only involves the value a of the field at the 
boundary in the diffusion term, it also involves a in the nonlinear advection 
(in f c ) and in the nonlinear stabilisation (in f b ). But further, it also involves 
time derivatives of the boundary value a. The reason is clearly that changes 
in a take time to advect and diffuse into the interior of the first few elements 
and so the effect of changes in the boundary value a upon the evolution of 
the grid values will lag, as seen by the opposite sign of the coefficients in 
the two columns of (14) . This lag increases with the element size h and 



hence accounts for the h? factor multiplying every d. Similarly, though not 
recorded above, higher order analysis shows each second derivative appears 
only in the combination h A a. Such effects are important when we try to use, 
in large scale problems, the expected accuracy and stability of these holistic 
discretisations on coarse grids. 



4 Apply Neumann boundary conditions at a 
midpoint 

Consider PDEs with Neumann boundary conditions of prescribed gradient of 
the field u. Two different sorts of numerical approximations are commonly 
developed for such a boundary condition: a grid point is placed at the bound- 
ary (as for the Dirichlet problem of §|3|); or the boundary is arranged to be 
midway between two grid points. For the first case I found that the dis- 
cretisation of the diffusion is expressed in terms of non-symmetric matrices. 
For the second case, the matrices are symmetric which again corresponds to 
the self-adjoint nature of the diffusion operator. Thus I report here on the 
second case where the boundary is at a midpoint of the grid. 

Without loss of generality, let the Neumann boundary condition of pre- 
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scribed gradient at a grid midpoint beQ 

du 



h 



dx 



ja(t) at x = X1/2 



(17) 



where, as for the ibcs, we actually are interested in the case 7 = 1. To 
supplement this left-hand boundary condition on the leftmost element we 
use the ibc (f|) as before. We execute the computer algebra program of Ap- 
pendix |A| with options dirichlet : =0 and midpoint : =1 .0 Again the interior 
discretisation is whereas in elements the near the boundary we find the 
grid values evolve according to the form (|T0|) where now the matrices are 
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and G3 is as in (|13D . Denote by D the symmetric matrix appearing above in 
T> linear in 7 and then row extended across the interior of the domain: D is 

Neumann conditions at a grid point, say X\, may be incorporated into the analysis by 
similarly requiring hdu/dx = 7a at x = x\ in place of the left-hand ibc (||). 

2 The computer algebra algorithm takes a few more iterations to complete because 
the algorithm is tuned to the discrete form of the ibc (H) whose left-hand side is only 
approximately that of (fTi 
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the matrix of the second-order centred approximation to the second deriva- 
tive, S 2 , but incorporating a zero derivative condition to the left of the grid. 
Observe in the order 7 2 and 7 s matrices of T> that the discretisation near the 
boundary of the diffusion term, fllTf ), is simply the interior discretisation seen 
in the first line of (f|) with the truncated D replacing the centred difference S 2 . 
Although here there is no clear pattern in C nor £>, as written above see that: 
B is numerically close to D/12; and, upon denoting C as the matrix linear 
in 7 in C, the 7 2 matrix is (DC + CD)/ 2 just as for the Dirichlet boundary 
conditions. However, in this case the identification is a little forced because 
numerically small discrepancies have been gathered into g c along with in- 
distinguishable similar terms corresponding to -^(5 2 Uj n5 3 Uj + 5 4 Uj fxSuj) 
appearing in the interior discretisation (||). Nonetheless the closeness of the 
match and the appropriate symmetry is reassuring. 

The forcing from the boundary takes the form 
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As discussed in the previous section, the boundary value for the flux a appears 
in a wide range of terms in the discretisation near the boundary. The reason 
again is that the flux feeds into the subgrid scale fields of the boundary 
elements and so affects the interaction between the various physical processes. 
The scope for such interaction increases with increasing element size h and so 
accounts for the appearance of the h 2 factor in front of the time derivative a. 
We need to know such effects on coarse grids. 
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5 Conclusion 

We have considered the discretisation of Burgers' equation ([J) on a finite 
domain as a worked example of incorporating physical boundary conditions 
into the derivation of holistic discretisations. The two most common physi- 
cal boundary conditions were considered in §[3] and §£|. Crucially we found: 
it is easy to maintain the symmetry appropriate to self-adjoint differential 
operators; discretisations are developed with constant bandwidth across the 
whole domain; and lastly our resolution of subgrid scale processes shows 
how time derivatives of the boundary forcing should also be included in the 
discretisation. 

Additionally, using the discrete form of the inter-element boundary con- 
ditions we observe for the first time a high order consistency of the 
nonlinear dynamics of Burgers' equation. 
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A Computer algebra derives different bound- 
ary discretisations 

Straightforward computer algebra programs are written to find the centre 
manifold and the evolution thereon 0, e.g.]. To ensure correctness and to 
provide a basis for further work I include the computer algebra code. This 
replaces extensive recording of elementary algebraic steps in the derivation 
of the results. 

I implement the construction algorithm in REDUCE^ The overall plan of 
the algorithm is to iteratively satisfy Burgers' equation ([1]) and the relevant 
internal and actual boundary conditions. A general interior element 

is analysed in lines 47-53 while the o : =3 elements near the boundary are 
analysed in the loop of lines 54-71. Although there are many details in the 
program, the correctness of the results are only determined by driving to 
zero (lines 53, 70 and 73) the residuals of Burgers' equation in each element 
and the internal and domain boundary conditions: lines 47-9 evaluate the 
residuals for an arbitrary interior element; lines 55-8 for near domain bound- 
ary elements; and lines 60 or 62 the domain boundary condition. The other 
details, such as the updates computed in lines 50-2 and 67-9, only affect the 
rate of convergence to the ultimate answer. 



1 Comment Find the discretisation of Burgers' equation, or 

2 variants, with three different boundary conditions: 

3 * Dirichlet, u=a(t) at x=x_0 

4 * Neumann, midpoint, u_x=a(t)/h at x=x_{l/2} 

5 * Neumann, gridpoint, u_x=a(t)/h at x=x_l 

6 (c) Tony Roberts, May 2001; 

7 % set options to 1 (true) or (false) 

8 dirichlet : =1 ; midpoint :=0; 

9 o:=3; % is number of near boundary elements 



3 At the time of writing, information about redu ce was available from An thony 
C. Hcarn. RAND, Santa Monica, CA 90407-2138, USA. |mailto : reduceQrand . org] 
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11 % improve printing 

12 on div; off allfac; on revpri; factor gam,h; 
13 

14 % make function of xi=(x-x_j)/h 

15 depend xi,x; 

16 let df (xi,x)=>l/h; 
17 

18 % solvability condition 

19 operator solg; linear solg; 

20 let { solg(xi-~p,xi)=>(l+(-l)-p)/(p+2)/(p+l) 

21 , solg(xi,xi)=>0, solg(l ,xi)=>l }; 
22 

23 % solves v"=RHS s.t. v(0)=0 and v(+l)=v(-l) 

24 operator solv; linear solv; 

25 let { solv(xi~~p,xi) => 

26 ( xi-(p+2)-(l-(-l)-p)*xi/2 )/(p+l)/(p+2) 

27 , solv(xi,xi) => (xi~3-xi)/6 

28 , solv(l.xi) => (xi~2)/2 >; 
29 

30 % parametrise with evolving u(j) and boundary forcing a 

31 operator u; depend u,t; 

32 let { df (u(~k) ,t)=>sub(j=k,gj) 

33 when (not f ixp(k) ) or (f ixp(k) and k>o) 

34 , df (u(~k) ,t)=>g(k) when fixp(k) and k<=o and k>0 >; 

35 operator a; depend a,t; 
36 

37 % linear solution in jth element 

38 array g(o) ,v(o) ; 

39 vj:=u(j); 

40 gj:=0; 

41 for j:=l:o do v(j):=u(j); 
42 

43 % iterative refinement to specified error 

44 % here work to error |u|~4+gam~4 by multiplying advection by gam 

45 let { gam~4=>0, df(a,t,2)=>0 >; 

46 repeat begin % first do interior elements 

47 deq:=-df (vj , t) -gam*vj *df (vj ,x)+df (vj ,x,2) ; 

48 rbc : =- (sub (xi=+l , vj ) -sub (xi=0 , vj ) ) +gam* (u( j + 1) -u( j ) ) ; 

49 lbc : =-(sub(xi=0 , vj )-sub(xi=-l , vj ) )+gam* (u( j ) -u( j-1) ) ; 

50 gd:=(rbc-lbc)/h"2+solg(deq,xi) ; 

51 gj:=gj+gd; 

52 vj :=vj+h~2*solv(-deq+gd,xi)+xi*(rbc+lbc)/2; 
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53 ok:=if (deq=0) and(rbc=0) and(lbc=0) then 1 else 0; 

54 for j:=l:o do begin % near boundary elements 

55 deq:=-df (v(j) ,t) -gam*v(j ) *df (v(j) ,x)+df (v(j) ,x,2) ; 

56 rbc:=-(sub(xi=+l,v(j))-sub(xi=0,v(j)))+gam*(u(j+l)-u(j)) ; 

57 if j>l then % internal left BC 

58 lbc:=-sub(xi=0,v(j) )+sub(xi=-l,v(j) )+gam*(u(j)-u(j-l)) 

59 else if dirichlet then % dirichlet at x_0 

60 lbc : =-sub (xi=0 , v ( 1 ) ) +sub (xi=-l , v ( 1 ) ) +gam* (u ( 1) -a) 

61 else if midpoint then °/ neumann at x_{l/2} 

62 lbc:=-sub(xi=-l/2,h*df (v(l) ,x))+gam*a 

63 else begin % neumann at x_l 

64 abc:=-sub(xi=0,h*df (v(l) ,x))+gam*a; 

65 lbc : =2*abc-rbc ; 

66 end; 

67 gd:=(rbc-lbc)/h"2+solg(deq,xi) ; 

68 g(j) :=g(j)+gd; 

69 v(j) :=v(j)+h"2*solv(-deq+gd,xi)+xi*(rbc+lbc)/2; 

70 ok:=if (deq=0)and(rbc=0)and(lbc=0) then ok else 0; 

71 end; 

72 showtime; 



73 end until ok=l; 
74 

75 end; 



AJ Roberts, February 1, 2008 



